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Nuclear thermal to electric power conversion carries the promise of longer duration 
missions and higher scientific data transmission rates back to Earth for a range of missions, 
including both Mars rovers and deep space missions. A free-piston Stirling convertor is a 
candidate technology that is considered an efficient and reliable power conversion device for 
such purposes. While already very efficient, it is believed that better Stirling engines can 
be developed if the losses inherent in current designs could be better understood. However, 
they are difficult to instrument and so efforts are underway to simulate a complete Stirling 
engine numerically. This has only recently been attempted and a review of the methods 
leading up to and including such computational analysis is presented. And finally it is 
proposed that the quality and depth of Stirling loss understanding may be improved by 
utilizing the higher fidelity and efficiency of recently developed numerical methods. One 
such method, the Ultra HI-FI technique is presented in detail. 

I. Introduction 

P OWER conversion with free-piston Stirling engines 1 - 2 has the potential to deliver high efficiency, low mass 
solutions for longer and more varied space missions. 3 - 4 NASA's goal of achieving 35% system efficiency 
(approximately 45% convertor efficiency) will require increasing the percentage of Carnot efficiency that 
Is achievable to 55%. In addition to using advanced high-temperature materials to increase the Carnot 
temperature ratio, it Is anticipated that advanced computational fluid dynamics (CFD) will help to identify 
the following losses ,,-,s (also shown in figure I): 

1 . Inefficient heat exchange and pressure loss in the regenerator, heater and cooler, 

2. Gas spring and working space loss due to hystcrisis and turbulence, 

3. Appendix gap lasses due to pumping and shuttle effects, 

4. Mixing gas losses from nonunifortn temperature and flow distribut ions perpendicular to primary engine 
flow axis, 

5. Conduction losses from the hot to cold regions 
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(>. Losses clue to combined radiat ion, conduction and convection in hollow spaces 

7. And in general, inaccurate lass representations due to use of 1-D How design codes to account for 
flow and heat transfer through area changes (between components) where phenomena such as flow 
separation and jetting from tube into regenerator may occur. 

In addition, the following artificial numerical losses must be considered when computational simulations 
are performed (also shown in figure 1): 

1. Moving/deforming mesh losses from repeated low order flow field interpolations. 

2. Transient/Unsteady heat t ransfer and flow loss from inconsistent and inaccurate time discretization, 

Flow loss from low order approaches resulting in effectively adding artificial dissipation terms along 
sliding interfaces, at structured/unstructured grid interfaces, and within interior. 

Minimizing those artificial losses is best accomplished through higher order approaches. 9 While this 
approach is common in aeroaconstics, computational electromagnetics, and exterior flow problems, high 
order techniques have not yet been applied to simulating a Stirling device. Moreover, the following difficulties 
are often encountered when using high-order approaches: 

1. Generation of high-order, smooth, body-fitted grids around complex configurations can be difficult. 10 

2. High-order formulations can lack nonlinear robustness. 10 

.'{. The general usefulness of high-order methods is limited by fiist order accurate shock capturing. 11 

Fortunately, with the exception of the possibly random geometry in the regenerator, the free-piston 
design is essent ially smooth and admits curvilinear structured grids. The issue of nonlinear robustness is an 
open Issue, but a proposed approach is detailed later. And finally, the working gas is subsonic and shockless 
throughout the entire region. 12 For these reasons, a high-order approach should be investigated for "whole 
engine" Stirling analysis. 

This paper will highlight the past and present Stirling analysis techniques and discuss their relative 
merits. 
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Schematic with Springs and Dampers 



Working Fluid (bounce space and w/orwng space) 


Figure 1. Schematic of Actual and Artificial Numerical Losses In a FYee- Piston Convertor 


II. Description of the Problem 

The dual opposed configuration sliowu in figure 2 I3.i7.is is being developed for nmltiinission use <i.e\. for 
use in atmospheres and space), including providing electric power for potential missions such as unmanned 
Mars rovers and deep space missions. 19 



Figure 2. Dual Oppused Stirling Convertors Reduce Vibration (Schrelber) 

An example of the desired full Siirliug engine simulation Is shown in Fig. 3. 16 The pseudo desigu Is 
presented to avoid export control issues. It is anticipated that full 3D simulations will provide a level of 
geometric and flow detail necessary for further design improvements. For example, hardware experiments 
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lmvo shown that largo performance gains can be made bv varying manifolds and heat exchanger designs to 
improve flow distributions in the heat exchangers. 13 ’ I4, Ir * 

The kinds of How and geometry that occur in a Stirling engine for 
which a modeling technique must address are as follows: 8, 1 

1. Oscillating flow which changes the effective conduction, flow frict ion 
coefficients and heat diffusion. 211 

2. Low much number How (no shocks). 

Compressible flow due to enclosed varying volumes and heat transfer 
effects, 

4. Laminar. Transitional, and Turbulent How with Reynolds numbers 
from UK) to 1001)0 (based on various length scales pertinent to the 
flow region), 

5. Conjugate heat transfer, thermal dispersion and local thermal 
nonequilibrium in the porous media regenerator, 

(>. Micron to Millimeter scale geometry 21 

7. Sliding Interfaces in the appendix gap 

8. Deforming Flow Regions Due to Compression and Expansion 

III. Stirling Cycle Analysis Approaches 

The techniques of analysis fur Stirling engines can be categorized with Martini’s 22,23 nomenclature as 
follows: 

A. Zeroth Order Analysis 

William Beale observed 21 that most modern engines operate under similar conditions of the parametric 
ratios: dead volume ratio, temperature ratio, swept volume ratio, displacer- piston phase migle advance. 
Thus correlating the experimental data available on various engines, he fotmd that one could represent the 
engine performance as: 



Figure 3. Example Transient 
Simulation of Pseudo Stirling 
Convertor-Colored by Velocity 
Magnitude 


P 0 = UMSp ntcan V, w f 


( 1 ) 


B. First Order Analysis 

The first analysis of Sliding engines was done in 1871 by Gustav Schmidt in which he obtained closed-form 
solutions for the special case of sinusoidal volume variations mid isothermal hot and cold spaces. 23 In Fig. B 
(a), the model equations from the isothermal assumption shown in Fig. 4 is shown. These can be integrated 
with a computer but a closed form solution Is available in which sinusoidal motion Is assumed. 

The relat ive volume changes can be calculated using the phasor diagram and schematic in Figs. 5 (b and 
c). This analysis is loss-free and a correction factor is normally applied in practical designs to account for 
the losses. 
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Figure 4. Ideal Isothermal Model (Uriell) 
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(a) Ideal Isothermal model set of equations 


(h) Volume Vector Diagram 


Figure 5. Closed Form Solution Process (Uriel!) 
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C. Second Order Analysis 

This level of analysis may be based on an adiabatic 
analysis that subtracts losses caused by beat trans- 
fer auil flow power losses. It. relics on a modified 
Schmidt analysis (Schmidt analysis normally implies 
isothermal hot and cold spaces) .uni requires nonlin- 
ear time Integration of the model equations. It as- 
sumes adiabatic expansion and compression regions 
as shown in Fig. 0(a). 

I. Heal Transfer Losses 

Heat transfer loss occurs throughout the system but 
notably in the appendix gap region (space between 
displacer and heat excliaugers) through the follow- 
ing processes:- 1 

• Shuttle heat transfer in which a conduction 
loss down the displacer or piston and associ- 
ated cylinder is enhanced by one wall's oscil- 
latory motion"' 

• Gas enthalpy transfer (pumping) which is the 
net enthalpy transport down the gap by virtue 
of the working gas motion, pressure and tem- 
perature. 

• Hysteresis heat transfer which Is the net heat 
transfer from the gas to the cylinder due to 
temperature gradients set up by the pressure 
variations in t he expansion space. 

In addition, conduction losses through the wall 
aud displacer must he accounted fur. 

H. Flow Power Losses 

Flow power loss occurs mostly in the regenerator re- 
gion due to pressure drop across the matrix (which 
may be enhanced by non-uniform flow). The vis- 
cosity of the working gas mid gas spring hysteresis 
losses must also he accounted for. 

,V. Berchawitz Dynamic ami Thennondynumic 
Closed Form 

A closed form expression (Fig. 7( which included 
masses for dynamic as well as thermodynamic anal- 
ysis was derived by Berchowitz. 2 " 



{.i I Ideal Adiabatic Model 


p = M R ! 0* /TV + Vk/T*-*- Vt + Vh/TTi-t- Vr ,'TV) 

-vp(dVc;Tck -d Verne) Pic«inc 

• [Vt ! Tck - , <Vk JTk * Vi m * Vh Hb) 

-Vc/Thc] 

n ii n ii n 

•d-a-s-a-o 

Mfssrs 

dmc _ .p dVt + Vc dp / r ) ) (P Tck) 

dmf : lip d"r - 'h dp / y) ) (R TV) 
drak - rak dp / p 
drai-uudp/p 
draii - raiidp / p 

Mass 

Accuinuiatiuis 

:rcK - -dne 
mkr - nick 1 - (bn:-: 
into - dice 
nuC’ - nuir i dmh 

Mass Fbw 

It .nick' > j Him TVk - Tc tls* TCk - T): 
if ml*' > Q dun Tie - Tli d* TV* - TC 

Cnrch'ional 

TViupeiatiues 

dTc - Tc (lp 1 p + dVc t Vt - «liuc ; ui) 
dTe - TV (dp ) p * dVc t Vc iliuc ) an) 

TVmpnaUues 

dQk - Vt dp cv / R - cp fTck lixk Tk ntta*) 
dQr - Vi dp cv/R- cp’ (7s uflo* - Tli null) 
dQh - Vli dp cv / R - cd (Huiali Tto rale') 
dWc - p «1V: 
dWe - p dVe 
dW- d\Vc + dWe 

W-lirWl 

Energy 


1 1>) Ideal Adiabatic Model Equations 


Figure 6. Adiabatic Solution Process (Uriel!) 
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D. Third Order Analysis 

Third order analysis uses control volumes or nodes 
to directly solve one-dimensional governing equa- 
tions. Some of the first analysis at this level 
of fidelity was by Finkclstoin , 28 Urieli , 29 and 
Berchowitz . 30 

The codes by David Gedeon referred to 
as G LIMPS ' 1 - 1,2 and Sage®- 34 -* are one- 
dimcnsional and solve the governing equations 
implicitly in space and time. The grid in- 
cludes all time because a periodic solution is as- 
sumed/forced. Therefore, it is not possible to 
model transient startup Ireliavior. 

The linearized harmonic analysis code re- 
ferred to as HFAST 30 solved a steady-state pe- 
riodic problem in the frequency domain. Again, 
transient behavior is not modeled. 

The code by Mart ini engineering 1, was never 
validated but claimed transient modeling capa- 
bility. It. is not clear the governing equations are 
being solved since it appears many simplifications 
based on experimental correlations are used. 

Another unvalidated but interesting code by 
Renfroe 38 attempted a one-dimensional analysis 
using explicit Runge-Kutta time stepping and a 
Newton solver to solve the nonlinear equations. 

Finally, the Stirling dynamic model (SDM ), 39 uses a one-dimensional analogy of an entire Stirling con- 
vertor by linking together representative elements within the Simplorer(TM) commercial software package 
by Ansoft Corp. This tool enables approximate whole convertor dynamic analysis. Recent work attempts 
to incorporate thermodynamics via David Gedeon ’s Sage code described earlier. 

E. Fourth Order (Multi-Dimensional) Analysis 

At this level of analysis, relat ively little has been completed because t he third order analysis is faster and for 
the most part lias been an adequate engineering tool. However, to improve efficiency further (to understand 
and reduce losses) it will likely require a better understanding of the actual flow features and heat transfer 
throughout the engine. 

/. Modified Computer Aided Simulation of Turbulence (CAST) 

The modified CAST code 40 is based upon the Semi-Implicit Method for Pressure-Linked Equations SIM- 
PLE 41 method but is restricted to two dimensions . 42 It was modified to include oscillatory boundary condi- 
tions and conjugate heat t ransfer. It has been used to model Stirling components but has not been extended 
to a whole engine simulation tool. 

A pressure-splitting technique was added 43 to reduce the computational requirements. It was based on 
separating the thermodynamic and hydrodynamic pressures so that these widely varying scales could l>e 
solved with less round-off error and better efficiency. 


Free-Piston Stirling Dynamic Equations 
(B*!Cl*>*<1U, 1980*) 
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Figure 7. Closed Form Adiabatic Solution (Schreiber) 
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2. CFD-ACE 


This commercial code has been used to model a two-dimensional representative Stirling engine. 44 45 ' It is 
also based upon the SIMPLE technique. The regenerator is not currently modeled correctly since thermal 
equilibrium Ls assumed between the gas and solid. The code does not currently support sliding interfaces on 
parallel computers which is critical for modeling the appendix gap region. 

This finite volume code can utilize both structured and unstructured grids. 

$. Fluent 

This commercial code is also based upon the SIMPLE method (and PISO method for high speed flows). It 
currently has similar regenerator modeling limitations in that it is designed for non-oscillating flows. It docs 
however have a sliding interface that could be used for appendix gap modeling on parallel computers. It 
Ls l>eing used by several commercial manufacturers for this purpose (references are proprietary). It is also 
finite volume based and can utilize both structured and unstructured grids. 

4 . STAR-CD 

The Simulation oflhrbulcnt Flow in Arbitrary Regions (STAR) 46 code also uses SIMPLE and PISO methods. 
Its companion product (STAR-HPC) Ls the parallel computer version. It also has sliding interfaces and 
deforming mesh capability. It has been used in the related field of internal combustion engine piston modeling, 
and some Stirling engines have been modeled with it. 4, 

5. CFX 

This code also used SIMPLE and PISO methods on unstructured grids. It also has sliding interfaces imple- 
mented, but no Stirling engine modeling with this software has been publicly published. 48 

6 . Others 

While there are other in-house codes, they are usually limited to modeling only specific regions of the Stirling 
engine such as the regenerator, or the displacer. 

IV. Whole Engine Modeling 

The use of two-dimensional C'FD models cart significantly extend the capabilities, compared to third-order 
analysis, for the more detailed analysis of the complex heat transfer and gas dynamical processes which occur 
in the internal gas circuit/’ 2 A standard k — c turbulence model has been used for such two-dimensional 
simulations, but the flow is constantly transitioning from laminar to turbulent and back. A better numerical 
model would be large eddy simulation due to the transient nature of the flow. More recently, full 3-D 
calculations have been performed with a commercial code. 53 The temperature results are similar to the 
second order method results. 5,4 The multidimensionally computed power, however, was about half of the 
second order prediction. Moreover, along the axis of the compression space it was found the change of the 
temperature of the working gas was quite different from harmonic in time. Third order modeling efforts can 
resolve higher harmonics of periodic time functions and therefore if enough grid points are available in time 
they can achieve similar results. However, turbulence is inherently a-periodic and it is not clear if third order 
modeling efforts which assume periodic solutions can fully capture this. 
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V. Regenerator Modeling 



A very important specific area of modeling difficulty is the 
regenerator. As shown in Fig. 8. since the regenerator (depend- 
ing upon one's definition of effective) lias roughly 3 to 40 time* 
more effective heat transfer than the heater/ 15 any inefficiency 
of the regenerator represents a significant loss for the entire 
Stirling engine. Hence, any numerical losses in this region will 
disproportionately influence the entire Stirling simulation. A 
new low loss numerical method will he proposed later to ad- 
dress tliis. 

A. Manifest 

Tlie code referred to as Manifest 31 solved the porous medium 
model equations shown below: 


Figure 8. Importance of Properly Modeling 
the Regenerator for System Studies (Uriel!) 


,/ , J O ( C* \ _ „ 

77T 1 TTZHTT \~r) ~ U 


-i(T.-w‘+5rf)-e=o 


&(xt.) + q = u 


( 2 ) 

where M = p.G = 0pV,E = pc, using the Beam and Warm- 
ing 50 implicit time stepping approach. This code was written 
in curvilinear coordinates and seemed to perform well for low Reynolds number jots impinging the matrix. 
It was not developed further because of excessive computational time. 




(c) Photographed Geometry 


Figure- 9. Regenerator Geometry 

Another important issue is the geometrical sliapc of the matrix in the regenerator. As shown in Fig. !), 
most regenerator models don't assume a precise geometrical shape for the elements of the regenerator. 
However, as reported in Park, 19 a simplc-to-fabricatc woven mesh, consisting of bonded laminates of two- 
dimensional plain-weave conductive screens can be manufactured to have a wide range of porosity and a 
highly anisotropic thermal conductivity vector. In addition to providing superior performance in many 


y 
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cases/’" the regular geometry greatly simplifies the analysis. However, those mesh sheets are not woven wire 
screens formed by etching used in regenerators which may in practice not be as precisely formed as mesh 
sheets. 


VI. Numerical Error and Efficiency 


The critical Issues in simulating a Stirling engine are: 

1. The artificial numerical losses should he a small percentage of the actual losses 

2. The simulation must complete in a reasonable amount of time. 


The first issue is a problem with many commercial codes that utilize low order methods. For example, 
this problem can be observed with a simple fiist order approximation of the first derivative: 


dF(t) __ AF _ F(t + At) - F(t) 
— ~£t At 


We may assume we are computing the derivative of a simple harmonic function. F(t) = t"~". 
the first order derivative approximation in Eq. (3) results in the solution: 


We find 


AF _ r. 
37" 


~X\ " 




lulr 


l-C( 






di 




Dissipation 


Dispersion 


("»> 


Higher order approximations of the derivative int roduce similar errors, but their relative effect is less. 


A. Time Advance Error 

In addition to the spatial and temporal discretizations which produce the dispersive and dissipative errors 
as discussed, another form of error which Is important for oscillating flows Is improperly modeled multi- 
dimensionsal boundary conditions. 

As an example, consider the pressure near the displacer wall. The numerical approach can be expanded 
as a Taylor series in time: 

p(r,„, ! + AI)= p(x,y, + + (5) 


If the time derivatives are written as space derivatives (for simplicity consider 2D Euler flows): 


dt 



dp 


du- 1 dy)) 


( 6 ) 


Then clearly these spatial derivatives must Ik? properly defined at a surface for the time derivative to be 
correct. Unfortunately, most commercial codes only impose the no-fiow condition and these spatial derivative 
terms are not properly posed at the boundary. ’ 1 This error represents the wrong time-accurate physics near 
walls. If the walls are moving (such as in a displacer) this error is more pronounced. 
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Higher time derivatives have many more terms in them that 
are generally incorrectly modeled ai surfaces. Fur example 
in Fig. 10. the colored terms represent commonly improperly 
posed terms. The pink represents first order or dispersive error, 
the green represents diffusive error, and the blue represents an 
error in the entropy terms. In au oscillating How loss study, 
these terms can become significant as time accumulates. 

B. Description of the SIMPLE method commonly 
used 

The Semi-Implicit Method for Pressure-Linked Equations 
(SIMPLE) is used in most of the commercial codes today. An 
understanding of its effect on Stirling .simulations is important 
for interpreting the losses found from present computational 
analysis. 



Figure 10. Second Time Derivative 


C. Improving Efficiency 

One of the primary limit at ions of using commercial codes today 
for three-dimensional analysis of Stirling cycle engines is the 

time it takes for second order accurate analysts. A recent trend wliick improves this situation Is relatively 
low-cost parallel computers (pie lined in Fig. 11(a)). The trend is for more processors on a single hoard and 
for liigher speed communications linking many boards together. Many commercial coties using the SIMPLE 
method can now utilize parallel computing, however, there still secuis to he a practical limit due to the 
implicit nature of the method. Further improvement in performance can be expected by utilizing high order 
methods as shown in Fig. 11(h). The higher order methods in coujuctiou with parallel computing offer the 
promise of design speed analysis. An explicit method with these properties Is presented in the next section. 



( 11 ) Cluster Computing 


Relative Tine Coot or Wave Resolution by Irrar ttequiraMOC 



|li) Efficiency Chart 


Figure 11. Enhancing Simulation Performance 
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VII. Proposed Computational Stirling Analysis Numerical Approach 


Computational Stirling Analysis (CSA) requires efficient, 
high-resolution simulation tools. Mast current CSA techniques 
utilize low order finite- volume approaches because high order 
accuracy is considered too difficult or expensive to achieve and 
hopefully unnecessary for some problems. As shown in Fig. 12, 
these low order techniques require a relatively dense grid com- 
pared to the high order approaches. Two notable exceptions 
are the finite-volume WENO and finite element Discontinues 
Galerkin Methods (DGM). 57 However, they tend to be inef- 
ficient computationally due to their use of unstructured grids 
and complicated mathematics. 

Although finite difference techniques have been developed 
which maintain low losses over time, they have great difficulty 
insuring flow conservation in crucial locations (i.e. appendix 
gap. moving geometry). Moreover, the finite difference form 
produces grid singularities at block corners of a multiblock 
curvilinear grid. In three dimensions, grid singularities become linos of only first order accuracy which 
swirling turbulent flow will force eonvcctcd disturbances through. 

It's accepted common practice to use low-order finite volume techniques in computational fluid dynamics 58 
precisely because of the drawbacks of finite-difference techniques in complex geometry. CSA has focused on 
low order finite-volume techniques because it is not currently clear how to extend finite-volume techniques 
to very high order. Only recently have modest improvements in the accuracy been achieved. For example, 
Agarwal 59 lias developed a fourth order finite volume technique for the Euler equations on unstructured 
grids. And, Wang 60 developed a fourth order optimized structured/curvilinear finite-volume technique. And, 
Rezgui 61 has developed a third-order finite-volume Euler solver. A low dispersion finite volume technique 
by Nance 02 used Tam's dispersion relationship preservation on the interpolated flux to improve spatial 
resolution. And most recently. Z.J. Wang 63 developed a spectral volume method in which more degrees of 
freedom are made available by using subcells within a cell. As shown in Figs. 13 and 14, the third order 
method produced better results while using exactly the same number of degrees of freedom. 

These recent developments demonstrate the feasibility and importance of high order low dissipation low 
dispersion finite-volume techniques. But their basic approaches severely limit their potential for very high- 
order extensions which arc known to produce significant performance improvements' 51 without the need for 
special opt imizations. 

Building upon previously established Hermitian finite-difference approaches, 51 this work utilizes the spa- 
tial derivatives of the convective fluxes to insure exceptionally high resolution and efficiency in complex 
geometries. The technique can be used to model a complete Stirling engine. 

VIII. Numerical Approach 

For pedagogical simplicity, the two-dimensional Euler equations will be shown, but the extension to fully 
viscous flows Ls straight-forward. The finite- volume form of the Euler equations is normally solved in physical 



Figure 12. Points Per Wavelength Compar- 
ison 
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(b) 

(C) 

(a) 2 nd Older Faceted 


(b) 2 nd Older Curved 

(c) 3rd Order 

Figure 13. Density Contour-Comparing faceted geometry (a), curved geometry <b), and curved geomtry higher 
order (c) 
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Figure 14. Entropy Production-Comparing faceted geometry (a), curved geometry (b), and curved geomtry 
higher order <c) 
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However, recently structured multiblock grid generators have become available that can generate a very 
smooth curvilinear grid. Therefore, the following finite-volume form of the Euler equations can be used 
with many advantages over the traditional physical coordinate system finite volume (used in commercial 
packages) and the curvilinear multiblock finite difference methodologies . 67 - 68 


pv 

puV + n x p 
fivV + n y p 
pHV 


, V = v • H = n x u + riyi' 
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The jacobian determinant, ./, is defined by: 
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And the original Cartesian solution and flux vectors are: 
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( 11 ) 


E is total energy (internal energy plus kinetic energy) per unit mass, E = e + 

The surface integrals can be simply performed in the computational space since all cells have unit spacing. 
Since derivatives of the solution are also available one can assume a polynomial shape function 69 and the 
integrals become: 


fo 

fo'GifonW 


f («oo + a 10^ + O 111'/ + Snfyfidi) 

riooAi/ + flio^Ar/ + 2bij(Af}) 2 + «! , AC(Ar/) 2 

/ <«oo + 2io{ + ao\i) + «ii<//)r/< 

•TuuAi + a w i (AO 2 + «oi (AO'/ + anj (A*) 2 '/ 


( 12 ) 


where A? = Ai/ = 1. The fljj coefficients can be found using Hermit ian divided differences 51 As was done by 
Agarwal 59 and similar to Goodricli 70 the governing equations may lie differentiated as follows in curvilinear 
coordinates: 

/„(«&)•**» = -/«,(»<-»«) («> 
The fluxes may be calculated either using Jameson 7 1 style centered, but high order flux differences or by 
using a pseudo one-dimensional Riemann solution between the left and right states as follows: 59 


p. = I 

1 m.n 


+ Pirn,, -S|A|A-'S- 1 (/' lm ,, - F lm .„) 

L High! Left night Left 


(14) 


where S is the similarity matrix for diagonilizing the Jacobian <)F\ m „ft MF m-n , A is the eigenvalue diagonal 
matrix, and |A| refers to the spectral radius of A. 

With these developments many arbitrary accuracy finite volume schemes can be implemented to arbitrary 
order accuracy in space. Tin* time advancement may be implemented with Runge-Kutta schemes for low 
frequency Stirling simulations that contain high wavenumber requirements (steep flow gradients). Or the 
recursive time advance technique 51 could be used to achieve arbitrary accuracy in time as well. 


IX. Curvilinear Grid Issues 

One of the primary advantages of this finite- volume formulation over the common finite difference ap- 
proach is in resolving grill singularity issues as shown in Fig. 15. 

For the following reasons, a structured curvilinear multiblock arbitrary accuracy finite volume CSA code 
should be developed: 
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(a) Metric Term Gradient Distributiuu 


(b) Grid Singularity Avoided With Fimte Volume 


Figure 15. Grid singularities are unavuidable in complex geometries and act as a source of 1st order noise for 
finite difference methods. For example, calculating a derivative with a 7-point stencil (marked with "A" or 
”D") Is indeterminate. A finite volume (with volumes "V”) method communicates through faces and hence 
this grid poses no problems for It. 


• Calculating spatial derivatives is reduced to V order accuracy at grid singularities with a fiuite- 
difTeroiicc method. However, a finite-volume method communicutes through cell faces only and for the 
grill singularity shown all faces are tuiiquely defined. 

• In three-dimensions this advantage is significant because grid singularities become lines of first order 
error. The swirling turbulent How in the expansion and compression spaces will cause all couvoctcd 
disturbances to travel through this region with many undesirable side effects. 

• In addition to enabling high accuracy everywhere this formulation will resolve high curvature without 
requiring ns many grid points when the fluxes are written in strong conservation form. 

• Parallel implementation is direct as well and docs not require changing the wave propagation speeds 
that converting from an implicit compact to an explicit centered scheme dot*. 

• Flat plates are easily solved and therefore code validation can be performed on standard oscillating 
flow flat plate heat exchanger before trying more complex problems. 

X. Time-Stepping Options 

One significant difference with time-stepping finite- volume cells relative to finite-difference and finite- 
element cells Is the time derivatives arc actually cell averages. So it Is possible to contain physical shocks 
within a cell mid its cell average derivatives still are defined. Thanks to the property that integration ignores 
singularities, it seems plausible to have both high order accuracy and shocks without sacrificing either. 

Standard Ruiigp-Kutta time-stepping can lie used, but it’s not recommended because fidelity and stability 
is sacrificed since currently available methods effectively widen the stencil because they do not properly 
account for the presence of spatial derivative information on t he grid. 
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A. Arbitrary Time Accuracy 

If higher operating frequencies in Stirling devices are encountered then a recursive 51 technique of time- 
stepping can Ik' used on each cell. Again, since the mixed space-time derivatives are actually cell averages, 
they are st ill defined even in shocked flow. 

XI. Flux Interpolation 

It is necessary to convert the spatially-averaged cell quantities to a form for interpolating t he fluxes at the 
cell edges. The spatially averaged flow variables can he used as constraints as done by Lomax . ' 2 This requires 
solving a matrix system for a polynomial that when volume integrated matches the spatially-averaged cell 
quantity. If Hcrmite data is on tin* grid, then for example, a single cell interpolation with /, f,, f r , produces 
a fourth order technique if the left and right interpolated fluxes are averaged. However, by using data from 
two neighboring cells we can get b' h order accuracy. 

If we have up to s derivatives in each cell, as shown in the figure 10, 


j-1/2 

a a 

— Ar— 

1+1/2 

j-2 j-1 

L 

R ' L 

PH 


Figure 1G. Control volume in one dimension 

we can assume that the polynomial interpolating the flow variables looks like for cells j and j + 1 : 

2* + l 

“«> = £ “iC 

i=0 

The following constraints are applied to find the coefficients a,: 

i r-W 2 

37 J— Ax/2 ' W- "s - 

i r 3Ar / 2 

37 JAi/2 Hr"' "< TTP 7 ' 

with ni=0 up to s. 

The flux can then Ik- calculated at the cell faces using the high order interpolant. 

A 41*' order method using two cells would require inverting a 41 x 41 matrix which gets expensive. 
The problem is worse in multiple dimensions. Instead, we can use the extended divided differences 51 and a 
process known as reconstruction via the primitive function 64 or RP (see Colella and Woodward.™ 

Steps are for all m=() up to s: 

Step 1 Find F m (x _ ,, 2 ) = 0. F m (x l/2 ) = A xf mj F m (x 3/2 ) = A xf mj + A xf j+l 

Step 2 Use Hermitian divided differences on the data in Step 1 to create an interpolation polynomial. Choose 
a local coordinate system with £ = 0 somewhere in the two cells. The coefficients of P(£) = t\< + /<i£ + 
jf» 2 £ 2 + ... can be found. This function is also easily differentiated. 


( 15 ) 

(15) 
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Stop 3 And then, since the sample function f (£ ) = '.IL and F(x) = P(j‘), a function approximating f(x) is 

**> = * 3 * 

For constant Ar can use reconstruction via deconvolution. 

With t hese developments high order finite volume integrat ion is possible. The spat ial accuracy is high 
via the previous step. Now we can time advance by converting higher order cell averaged time derivatives 
to first order cell averaged time derivatives using the Cauchy- Kowalewski approach. 

So we time advance the averages (using the fluxes and their derivatives). We get new cell-averaged 
quantities at the next time level. We repeat the reconstruction to get new flux interpolants. 

XII. Conclusion 

A review of t he techniques used for analyzing the Stirling cycle lias been presented. .As the fidelity of this 
analysis has improved, so has the efficiency of the cycle in realized designs. Recent attempts at modeling 
a full Stirling engine have utilized commercial codes which have recently enhanced their ability to simulate 
moving, tieforming, and sliding meshes. Further improvements in the analysis will require addit ional research 
in the key areas of regenerator modeling and higher fidelity numerical approaches. 

Tin* regenerator is perhaps at the same time the most difficult and important region to model. A 
combination of direct numerical simulation, better theoretical representations, and high efficiency numerical 
methods on parallel computers will l>e required. Some of the more recent higher fidelity simulations of porous 
media were presented and appeared to indicate a trend towards modeling wit h higher order methods. 

Also, a review of the advantages of recently developed high fidelity finite volume methods was presented 
and a novel numerical approach for achieving extraordinary levels of fidelity in Stirling engine simulations 
was described. And, for completeness, the bibliography contains additional references that expand upon t he 
review presented but due to the constraints of space and time have not been specifically commented on here. 
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